%MacCormack scheme
function u=MacC(u0,dt,maxa,N,al,b)
u=zeros(1,N);
h=1/N;
la=dt/h;
flux = @(uu) uu.^al/b;
for j=1:N 
    if j~=1&&j~=N
        ustar=u0(j+1)-la*(flux(u0(j+1))-flux(u0(j)));
        fp=(flux(ustar)+flux(u0(j)))/2;
        
        ustar=u0(j)-la*(flux(u0(j))-flux(u0(j-1)));
        fn=(flux(ustar)+flux(u0(j-1)))/2;
    elseif j==1
        ustar=u0(j+1)-la*(flux(u0(j+1))-flux(u0(j)));
        fp=(flux(ustar)+flux(u0(j)))/2;
        
        ustar=u0(j)-la*(flux(u0(j))-flux(u0(N)));
        fn=(flux(ustar)+flux(u0(N)))/2;
    else
        ustar=u0(1)-la*(flux(u0(1))-flux(u0(j)));
        fp=(flux(ustar)+flux(u0(j)))/2;
        
        ustar=u0(j)-la*(flux(u0(j))-flux(u0(j-1)));
        fn=(flux(ustar)+flux(u0(j-1)))/2;
    end
    u(j)=u0(j)-la*(fp-fn);
end